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Abstract. We analyze the deviations from Maxwell-Boltzmann statistics found in recent experiments study- 
ing velocity distributions in two-dimensional granular gases driven into a non-equilibrium stationary state 
by a strong vertical vibration. We show that in its simplest version, the "stochastic thermostat" model 
of heated inelastic hard spheres, contrary to what has been hitherto stated, is incompatible with the ex- 
perimental data, although predicting a reminiscent high velocity stretched exponential behavior with an 
exponent 3/2. The experimental observations lead to refine a recently proposed random restitution co- 
efficient model. Very good agreement is then found with experimental velocity distributions within this 
framework, which appears self-consistent and further provides relevant probes to investigate the universal- 
ity of the velocity statistics. 

PACS. 45.70.-n Granular systems 05.20.GgClassical ensemble theory 51.10.+yKinetic and transport the- 
ory of gases 



1 Introduction 

Whereas equilibrium statistical mechanics has reached a 
rather mature phase, the understanding of non-equilibrium 
processes is far from complete. In particular, granular (and 
thus inelastic) gases £Q driven into a non-equilibrium steady 
state by a suitable injection of energy define a stimulat- 
ing research field where theoretical predictions can be con- 
fronted against model experiments, with the aim to under- 
stand the possible deviations from equilibrium behavior. 
A good probe to quantify these deviations is the velocity 
distribution of the grains, P{v), which has focused sus- 
tained attention recently, and has been shown to exhibit 
pronounced differences from Maxwell-Boltzmann statis- 
tics [2E1E1EK]- Several authors reported a stretched ex- 
ponential law [on the whole range of velocities available, 
which covers an accuracy of 5 to 6 orders of magnitude 
for P(v)} 

P(v)(xexp[-(v/vo) v ], (1) 

with an exponent v close to 3/2 3,4,6 (here v is the 
"thermal" r.m.s. velocity). This behavior was observed for 
the horizontal velocity components of a vertically vibrated 
2D system of steel beads in a wide range of driving fre- 
quencies and densities @], but also in a three dimensional 
electrostatically driven granular gas 0. 

At this point, some questions naturally arise, that will 
be addressed below: (i) is it possible to find a consistent 
model where this velocity distribution would emerge? (ii) 
what physical ingredients are required? 

One possible approach consists in performing "real- 
istic" molecular dynamics simulations. The model of in- 



elastic hard spheres (IHS) with binary momentum con- 
serving collisions, and a "reasonable" restitution coeffi- 
cient ITJ provides the simplest candidate. The energy loss 
in a collision is proportional to the inelasticity parameter 
1 — cvq where ao is the coefficient of normal restitution 
(0 < aa < 1), which in the simplest and efficient approxi- 
mation is a constant independent on the relative velocity 
of colliding partners. Such an approach has been presented 
in [715] ; and allows to reproduce the experimental veloc- 
ity statistics with a good accuracy. The possible lack of 
universality has also been addressed in [Hj. 

Another route, which contrary to the previous numeri- 
cal one has the merit to allow an analytical derivation of v 
in some cases JH], consists in formulating an effective mod- 
eling of the energy injection, considering idealized homo- 
geneous systems of inelastic hard spheres (given the exper- 
iments reported in 0], the assumption of homogeneity is 
well founded, see below). From this point of view, a simple 
and popular model consists in IHS with constant inelastic- 
ity, with a homogeneous forcing described by a "stochastic 

tracted attention, in particular because it has been shown 
analytically [5j that P(v) exhibits a high energy tail of the 
form of Eq. Q with v — 3/2, independent on dimension 
and restitution coefficient, in apparent agreement with the 
experiments. This result holds at the mean-field level for 
an homogeneous system. The above model, where an ex- 
ternal white-noise driving force acts on the particles and 
thus injects energy through random "kicks" between the 
collisions, is therefore considered to provide a relevant the- 
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oretical framework to quantify the non Gaussian character 
of velocity distributions. 

However, we shall see below that this uniformly heated 
model -in its simplest version- is unable to reproduce the 
experimental data: if simulated in dimension 2 or higher 
with an experimentally relevant value of the restitution 
coefficient (ao between 0.7 and 1), the obtained distri- 
bution P(v) is indistinguishable from a Gaussian within 
the experimental accuracy , in fact, the range of veloc- 
ities for which the high energy behavior exp[— (w/wo) 3 ^ 2 ] 
may be observed is even beyond reach of precise numer- 
ical procedures and corresponds to a regime where P(v) 
is practically vanishing [lower than 10 _6 P(0)]. The pre- 
dictions of this model are consequently incompatible with 
the experimental velocity distributions, that show impor- 
tant non Gaussian features already at thermal velocity 
scale jHElIEl- We emphasize however that generalizations 
of the aforementioned heated model have been proposed 
|17ll8j . With a convenient choice of the extra parameters 
introduced, one may obtain a velocity distribution close 
to that measured in the experiments (we will come back 
to this point in section |3J|. However, in such approaches, 
energy is injected in between the collisions whereas in the 
experiments we are interested in, the transfer of horizon- 
tal momentum takes place at every inter-particle collision 
(see section 0J . We will therefore focus on this feature 
for the 2D experiment reported in 0] and investigate in 
details the collision dynamics in the horizontal direction 
(section F2J). The vibrated system under study there shows 
important density and granular temperature gradients, es- 
pecially close to the boundaries which inject energy, but 
since the shaking is violent, there is a region where both 
gradients are very small simultaneously. The velocity ac- 
quisition in 0j has been restricted to this region, where 
the system, although open, may be considered as homo- 
geneous. In this article, we thus consider the following 
question : remaining at the level of a homogeneous sys- 
tem, what ingredients are required for a self-consistent 
effective description of the horizontal degrees of freedom, 
that exhibit the stretched exponential law JU ? 



2 Effective restitution coefficients 

In order to characterize the collision process in the hori- 
zontally projected system, we have measured directly the 
effective ID restitution coefficient from the experimental 
data provided by K. Feitosa and N. Menon, for a gas of 
stainless steel spheres (the system investigated in rj]) but 
also for glass, brass and aluminum beads JHI, which allow 
to sample a wide range of nominal inelasticities. Let us re- 
call briefly the experimental set-up. The balls (diameter: 
d = 1.600±0.002 mm) are confined to a vertical, rectangu- 
lar cage (32 d high x 48 d wide x 1.1 d thick) sandwiched 
between two parallel plates of Plexiglas. The cage is vi- 
brated vertically at a frequency of 60 Hz and amplitudes 

1 if this model is simulated in dimension 1 with ao £ [0.7; 1], 
the obtained non Gaussian behaviour at thermal velocities cor- 
responds to v > 2 instead of v ~ 3/2 in the experiments. 



up to 2.4 d, producing maximum accelerations, r, and 
velocities, Vq, of 56 g and 1.45 m/s respectively. The mo- 
tions of the balls are recorded with a high-speed camera 
which allows a location of each ball with a precision of 0.03 
d. The results we discuss here are taken in a rectangular 
(10 d x 20 d) window around the geometrical center of 
the cell, where, as mentioned above, density and granular 
temperature are almost homogeneous Moreover, the 
measured velocity distributions do not vary with height 
nor with the phase of the vibration cycle. The experimen- 
tal data can thus be considered as obtained in the bulk of 
a two-dimensional homogeneous (but open) system, rea- 
sonably far from the boundaries. 

The horizontal component of relative velocities are com- 
puted before (g x ) and after (g*) each collision, from which 
we deduce the effective restitution coefficient 



Fig- displays the histogram /x eX p(aid) obtained from the 
experimental data for different materials. At large aid, a 
power-law tail is evidenced. Note that values aid > 1 are 
expected, due to the transfer from vertical to horizontal 
translational kinetic energy FJOj • 

The strong correlations between relative horizontal ve- 
locity g x and aid are clearly seen in the scatter-plot (inset 
of Fig. rj) , with a very sharp cutoff above the second bi- 
sector aid oc l/gx- This cut-off follows from the definition 
of aid'- since the post-collision velocity is finite, large val- 
ues of aid may only be obtained for small values of g x , 
in which case Eq. (J2J) implies that the maximum aid is of 
order 1 / g x . 

These features are qualitatively the same for all mate- 
rials investigated. Moreover, three different densities have 
been investigated for steel beads, and the same distribu- 
tions have been obtained |19| . Similar distributions have 
also been measured in a different experimental set-up with 
rolling beads [2J ■ The relatively small number of collisions 
investigated does not however allow us to get accurate his- 
tograms for the joint distributions fj,(aid, g x ) nor the con- 
ditional fi(aid\g x )- The correlations evidenced in Fig. 
nevertheless play a crucial role, as will be shown below. 

More insight into the conditional distributions fJ.(aid\g x ) 
has been obtained in molecular dynamics of two-dimensional 
IHS driven by vibrating walls in Ref. 0. Histograms of 
fj>(aid) similar to the experimental results have been ob- 
tained. It turns out that the global /i(aid) is almost in- 
sensitive to the details of the system (density, velocity of 
the vibrating walls...), while this dependence exists for 
n(aid\gx) (and also for P(v)). The following character- 
istics of fi(aid\g x ) have been obtained: at constant g x , 
n(aid\g x ) is almost constant for aid € [0, aq], has a small 
peak at ao, and decreases as exp(— A{g x )a\ d ) for aid > 
a , with A(g x ) oc g 2 x . 

3 Random Restitution Coefficient model 

We have considered the possibility to mimic the experi- 
mental distributions by including randomness in the resti- 
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Fig. 1. Experimental histogram of aid for steel beads (circles) 
and glass beads (squares). Line: [i{a\d) obtained in the RRC 
model as selected by the collisional dynamics (see text). Inset: 
Experimental scatter-plot of aid versus relative precollisional 
horizontal velocity g x , for glass beads. Data from [1| and |19) . 



tution coefficient, following the approach of [201 where 
a random restitution coefficient (RRC) model was intro- 
duced to account for the fact that, in vertically shaken 
granular gases, the energy is transferred to the vertical 
degrees of freedom by the moving piston, and then to 
the horizontal ones through grain/grain collisions only. 
The heating of horizontal degrees of freedom thus occurs 
through the inter-particle collisions, and not in between 
as in the "stochastic thermostat" approach. Moreover, a 
globally dissipative collision may correspond to an en- 
ergy gain for the horizontal components of the velocities. 
This leads to the study of horizontally projected collisions 
with an effective restitution coefficient that can be either 
smaller or larger than 1, as the experimental data of Fig. ^ 
indeed show. In our case, the RRC model is therefore an 
effective approach in 1 dimension (since the original colli- 
sions are two-dimensional), in which IHS undergo binary, 
momentum-conserving collisions with a restitution coeffi- 
cient a±d drawn randomly at each collision from a given 
distribution [i(aid) which should mimic n e xp{ot-id) ■ Even if 
the original collisions are not random, the projected ones 
may be considered as such. 

Note that the energy injection is here given solely by 
the values of a\d larger than 1. The simulations are per- 
formed at the mean-field level of the homogeneous non- 
linear Boltzmann equation for point particles, solved by 
the Direct Simulation Monte Carlo method (DSMC) |22|. 
The velocity statistics is computed in the non-equilibrium 
steady state which is reached after a transient. 

We have first considered distributions with a large tail 
in order to reproduce the experimental n{ocid), but with- 
out any correlations with the relative velocities of the col- 
liding particles; in this case, it turns out that P(v) has a 
power-law decay at large v, in marked contrast with ex- 
perimental results 2 . This approach consequently needs to 
be refined and the next crucial step is to take into account 
the correlations between aid and g x , with the insight given 



2 It is noteworthy that a large v power law is actually ob- 
tained for any distribution /i(aid), see |20) . 



by the experimental scatter-plot (inset of Fig. Q and by 
the molecular dynamics results [7J- 

We have thus used distributions decaying as n(a\d\g x ) ~ 
exp[-(ai d g x ) 2 /R] at large au, where g x = g x /v is the 
rescaled velocity defined from the total kinetic energy of 
the system (v 2 = (v 2 )), and the parameter R can be varied 
with values of order 1. The function /J.(aid\g x ) is the only 
input needed to simulate the RRC model. For the consis- 
tency of the approach, the distribution n(aid) measured 
in the simulation needs to be close to its experimental 
counterpart. This comparison is displayed in Fig. ^ and 
justifies a posteriori the choice made for n(aid\g x )- Both 
experimental and numerical distributions /i(aid) display 
a power law tail of the form a^ with n ~ 3. 

The velocity distribution obtained from the RRC model 
is compared to the experimental measure in Fig. The 
agreement is satisfactory over the whole range of veloc- 
ities; in particular, the RRC distributions is compatible 
with the stretched exponential behavior reported in 0j, 
with an exponent v close to 1.5. Note that, since no pre- 
cise experimental data is available for the distributions 
of restitution parameters conditioned by relative precolli- 
sional velocity, the parameters R are tunable [as long as 
the global fj,(aid) coincides with the experimental one], 
and the value giving the best agreement has been chosen 
(2 < R < 4). The agreement remains satisfactory upon 
changing R, provided that the resulting large a cutoff re- 
mains sharp (i.e. R should not be too large). 

The RRC model therefore provides a self-consistent 
framework which allows to reproduce the experimental 
P(v) if implemented with the correct distribution of effec- 
tive coefficients. Moreover, the velocity statistics depends 
on the distribution of effective restitution coefficients: a 
broader [i(a\d\g x ) leads to a broader P(v), consistently 
with the numerical study of 7 which showed both broader 
P{v) and fJ-(aid\g x ) as e.g. the density is increased. Both 
distributions P and /x are equally sensitive to a possible 
non universality (dependence on material properties). As 
a consequence, an accurate experimental measure of \x 
appears as complementary to the direct computation of 
P(v), in order to assess the experimentally difficult ques- 
tion of the velocity statistics universality. 



4 Stochastic thermostat model 

For comparison, we have also considered heating of inelas- 
tic hard discs (2D) through the "stochastic thermostat" , in 
the framework of the non-linear homogeneous Boltzmann 
equation, where the large velocity tail has been shown to 
behave like exp(— v 3 / 2 ) this result has subsequently of- 
ten been considered as an agreement with the experimen- 
tal results. In this model, the energy injection is achieved 
through a random force rj (t) acting on each particle 

^=F + V (t), (r H (t)r )j (t'))=2D5(t-t% j (3) 

where D is the amplitude of the injected power and F the 
systematic force due to inelastic collisions. The variance 
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of r\ determines the granular temperature in the non equi- 
librium steady state, but has no influence on the form of 
the rescaled distribution function P(c x ). 

With the accuracy of Fig. [3] (the current experimental 
resolution), the corresponding numerical velocity distribu- 
tions are then found indistinguishable from a Gaussian, for 
physically relevant inelasticities in the range 0.7 < ao < 
1. Departure from Maxwell-Boltzmann behavior becomes 
manifest below ao = 0.6 (squares in Fig. EJ ^ which is un- 
physically low, but the velocity distribution is still incom- 
patible with its experimental counterpart. We also inves- 
tigated the possibility to describe the effective horizontal 
dynamics with the ID stochastic thermostat: for 0.7 < 
ao < 1 the velocity distributions are incompatible with 
the experimental P(v) displayed in Figs. [21 with opposite 
non Maxwellian features (underpopulated both at van- 
ishing and high energies More precisely, within the 
stochastic thermostat approach the kurtosis (v^) / (v 2 ,) 2 — 3 
of the distribution is negative for a > 1/V2 — 0.71 [5], 
irrespective of dimension, which corresponds to an un- 
derpopulated low velocity behavior at variance with the 
experimental data shown in Fig. [21 

We have also considered the stochastic thermostat mo- 
del for two-dimensional IHS with both tangential at and 
normal a n (— ao) restitution coefficients )24j . No numer- 
ical studies of P(v) can indeed be found in the litera- 
ture in this case, although an investigation into the non- 
equipartition between translational and rotational kinetic 
energies has been performed in |24| . The resulting veloc- 
ity distributions P{v) remain very close to a Gaussian 
for a n > 0.7 and arbitrary at (where —1 < a t < 1), 
as for smooth spheres (corresponding to a t = —1). How- 
ever, such a two parameter model may be too schematic 
compared to the experiments [21] and we have also con- 
sidered a more realistic approach with Coulomb friction 
along the simplifications discussed in [2Hj: a friction coef- 
ficient /i, is introduced in addition to (a ( ,a n ) |27j . This 
does not change significantly P(c x ) (see the squares and 
dashed line in Fig.^J. This seems to discard the relevance 
of such an approach for the comparison of the velocity 
distributions with experimental data. 

The above analysis shows that the stochastic thermo- 
stat in its original formulation (including some possible 
extensions) does not provide a relevant model of energy 
injection as far as the velocity distribution is concerned, 
although it may be useful to investigate other features 
such as kinetic energy non-equipartition in granular mix- 
tures |28|. However, a variant of this model may improve 
the picture. In particular, a multiplicative driving [corre- 
sponding to a velocity dependent amplitude D oc \v\ 2S in 
Eq. ©] has been studied in Refs. [T7HT%| . We have per- 
formed DSMC simulations for this model, choosing the 
value of S that, for a given reasonable inelasticity param- 
eter (a = 0.9), gives the best agreement with the experi- 
mental P(v). We obtained 6 ~ 0.6, which leads to the dis- 
tribution shown by the triangles in Fig. [21 The agreement 
with the experimental data is satisfactory, and displays 
a similar accuracy as obtained within the RRC model. 
Consequently, models with homogeneous energy injection 
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In conclusion, the Random Restitution Coefficient model 
with point particles captures the essential features respon- 
sible for the observed non Gaussian character of the ve- 
locity distribution P(v) in vibrated granular gases experi- 
ments pi], and represents therefore a self-consistent frame- 
work. The conditional distribution ^(aidl^) defining the 
appropriate collision rule has been shown to encode the 
relevant dynamic information and provides an alternative 
route to characterize the non equilibrium steady state, 
complementary to the direct measure of P(v). An inter- 
esting point would be to obtain an analytical prediction for 
fi(aid\g)- It is noteworthy that our approach is mean-field 
(Boltzmann) like, the only correlations considered being 
in the collision law. Its self-consistency, which was not an 
obvious point a priori, has been established by comparison 
with experiments. 

While we have restricted our analysis to the two-dimensional 
case, such investigations can be extended to three-dimensional 
systems, for which however experimental measures of ef- 
fective restitution coefficients seem more difficult. As im- 



Fig. 2. Rescaled distribution P(c x ) of horizontal velocities, 
on a linear-log scale (a) and linear scale (b). All distributions 
have the same variance (c^) = 1/2. Circles represent the ex- 
perimental data for steel beads 4': ; Filled triangles correspond 
to the Monte Carlo simulation of the RRC model, with a 1( i/gx 
correlations. 

may also describe quite accurately the experimental P(v), 
with the problem of predicting the values of the various 
parameters involved. 

5 Conclusion 
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Fig. 3. Rescaled distribution P(c x ) of horizontal velocities, 
on a linear- log scale (a) and linear scale (b). All distributions 
have the same variance (c x ) = 1/2. Circles represent the ex- 
perimental data for steel beads U; The squares correspond 
to a simulation of the "stochastic thermostat" with cto = 0.6, 
whereas for ao > 0.7, the corresponding P(c x ) is indistinguish- 
able from the Gaussian shown by the full line. The dashed line 
(very close to the squares) corresponds to a simulation of the 
three parameters model |25l26|2*7) with a n — 0.7, a t = 0.5, 
/i = 0.5. Filled triangles show the results for the multiplicative 
driving stochastic thermostat with 5 — 0.6 and a — 0.9 17 . 



plemented here, without an analytical knowledge of n(aid\g- 
the RRC model is not predictive since an experimental 
input is required to obtain the correct velocity statistics. 
Our results however suggest to assess experimentally the 
question of the universality of P(v) from the direct mea- 
sure of the distribution of effective restitution coefficients. 
These characteristics are indeed linked within the RRC 
model: the exponent and the whole shape of P(v), de- 
pend on the functional form of (J>(ttid\g x ). At this point, 
the fact that with a rather poor accuracy, similar P and /i 
have been obtained for the case investigated in |Hll9j sim- 
ply confirms the RRC picture, and calls for experiments 
with widely different collisional properties, such as hollow 
spheres. 
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